%function g=M(a0,b0,x0)

a0=1;
b0=0.1;
x0=80

% this code will not work if b0 is a negative integer
%if (b0==round(b0)) & (b0<0)
%	error('The code doesn't work if b is a negative integer');
%	return;
%end

% use identity:  M(a,b,x)=exp(x)*M(b-a,b,-x) 
if x0< 0
	x=-x0; a=b0-a0; b=b0;
else
	x=x0; a=a0; b=b0;
end

xmax=10;

% evaluate M(a,b,x) for x ge 0

%  use exact expressions for M if possible
if (x==0) | (a==0) 
	g=1
elseif a==b
	g=exp(x)
elseif (a==1) & (b==2)
	g=(exp(x)-1)/x

%  using the complete poly isn't the smartest idea if a is a very large 
%  negative integer but its done anyway 
elseif (a==round(a)) & (a<0)
	N=-a;
	g=1;
	cs=1;
	for in=1:N
		cs=cs*(a+in-1)/((b+in-1)*in);
		g=g+cs*x^in;
	end

%  otherwise, for small x use series definition
elseif x<xmax
	N=800;
	g=1;
	cs=1;
	for in=1:N
		cs=cs*(a+in-1)*x/((b+in-1)*in);
		g=g+cs;
		if abs(cs/g)<1.0e-12 
			break
		end 
	end

%  and use asym approx for large x
else
	N=40;
	g0=1;
	cs=1;
	for in=1:N
		cs=cs*(b-a+in-1)*(1-a+in-1)/(in*x);
		g0=g0+cs;
		if abs(cs/g0)<1.0e-12 
			break
		end 
	end
	q1=gamma(b)/gamma(a);
	g=x^(a-b)*exp(x)*q1*g0;
end

%  adjust answer if x0<0

if (x0<0)
	g=g*exp(x0);
end

g0
g




